############################################################################################################
# Allowed to Hold Public Office

m.m.d.1 <- lm(liberal ~ rel_banperm + rel_muslim + rel_dev + rel_rad + rel_immigrant + fem + I(age/10) + edu_high + religiosity + leftright, data=imp.data.1)
m.y.d.1 <- lm(rel_office ~ rel_banperm*liberal + rel_muslim + rel_dev + rel_rad + rel_immigrant + fem + I(age/10) + edu_high + religiosity + leftright, data=imp.data.1)

m.m.d.2 <- lm(liberal ~ rel_banperm + rel_muslim + rel_dev + rel_rad + rel_immigrant + fem + I(age/10) + edu_high + religiosity + leftright, data=imp.data.2)
m.y.d.2 <- lm(rel_office ~ rel_banperm*liberal + rel_muslim + rel_dev + rel_rad + rel_immigrant + fem + I(age/10) + edu_high + religiosity + leftright, data=imp.data.2)

m.m.d.3 <- lm(liberal ~ rel_banperm + rel_muslim + rel_dev + rel_rad + rel_immigrant + fem + I(age/10) + edu_high + religiosity + leftright, data=imp.data.3)
m.y.d.3 <- lm(rel_office ~ rel_banperm*liberal + rel_muslim + rel_dev + rel_rad + rel_immigrant + fem + I(age/10) + edu_high + religiosity + leftright, data=imp.data.3)

m.m.d.4 <- lm(liberal ~ rel_banperm + rel_muslim + rel_dev + rel_rad + rel_immigrant + fem + I(age/10) + edu_high + religiosity + leftright, data=imp.data.4)
m.y.d.4 <- lm(rel_office ~ rel_banperm*liberal + rel_muslim + rel_dev + rel_rad + rel_immigrant + fem + I(age/10) + edu_high + religiosity + leftright, data=imp.data.4)

m.m.d.5 <- lm(liberal ~ rel_banperm + rel_muslim + rel_dev + rel_rad + rel_immigrant + fem + I(age/10) + edu_high + religiosity + leftright, data=imp.data.5)
m.y.d.5 <- lm(rel_office ~ rel_banperm*liberal + rel_muslim + rel_dev + rel_rad + rel_immigrant + fem + I(age/10) + edu_high + religiosity + leftright, data=imp.data.5)


coefs <- cbind(coef(m.y.d.1), coef(m.y.d.2), coef(m.y.d.3), coef(m.y.d.4), coef(m.y.d.5))
c.combined <- apply(coefs, 1, mean)

var.within <- cbind(se.coef(m.y.d.1)^2, se.coef(m.y.d.2)^2, se.coef(m.y.d.3)^2, se.coef(m.y.d.4)^2, se.coef(m.y.d.5)^2)
var.within <- apply(var.within, 1, mean)

m <- 5
var.between <- (coefs - c.combined)^2
var.between <-  apply(var.between, 1, sum)
var.between <-  (m-1)^-1 * var.between

ses <- sqrt(var.within + (1 + m^-1) * var.between )

combined.res <- cbind(c.combined, ses)

combined.res.off <- rbind(combined.res, display(m.y.d.1)$n, display(m.y.d.1)$r.squared)
combined.res.off <- round(combined.res.off, 2)

colnames(combined.res.off) <- c("est", "se")
rownames(combined.res.off)[2] <- c("liberal_policy")
rownames(combined.res.off)[3] <- c("citizen_reaction")
rownames(combined.res.off)[13] <- c("liberal_policy*citizen_reaction")
rownames(combined.res.off)[14] <- c("N")
rownames(combined.res.off)[15] <- c("R2")

ACME.1 <- ACME.fun(m.m.d.1, m.y.d.1)
ACME.2 <- ACME.fun(m.m.d.2, m.y.d.2)
ACME.3 <- ACME.fun(m.m.d.3, m.y.d.3)
ACME.4 <- ACME.fun(m.m.d.4, m.y.d.4)
ACME.5 <- ACME.fun(m.m.d.5, m.y.d.5)

c.1 <- apply(ACME.1, 2, mean)
c.2 <- apply(ACME.2, 2, mean)
c.3 <- apply(ACME.3, 2, mean)
c.4 <- apply(ACME.4, 2, mean)
c.5 <- apply(ACME.5, 2, mean)

s.1 <- apply(ACME.1, 2, sd)
s.2 <- apply(ACME.2, 2, sd)
s.3 <- apply(ACME.3, 2, sd)
s.4 <- apply(ACME.4, 2, sd)
s.5 <- apply(ACME.5, 2, sd)

coefs <- cbind(c.1, c.2, c.3, c.4, c.5)
c.combined <- apply(coefs, 1, mean)

var.within <- cbind(s.1^2, s.2^2, s.3^2, s.4^2, s.5^2)
var.within <- apply(var.within, 1, mean)

m <- 5
var.between <- (coefs - c.combined)^2
var.between <-  apply(var.between, 1, sum)
var.between <-  (m-1)^-1 * var.between

s.combined <- sqrt(var.within + (1 + m^-1) * var.between )

combined.res <- cbind(c.combined, s.combined)
combined.res <- round(combined.res, 2)

av.ACME.off <- combined.res[3,1]
av.ACME.se.off <- combined.res[3,2]


ANDE.1 <- ANDE.fun(m.m.d.1, m.y.d.1)
ANDE.2 <- ANDE.fun(m.m.d.2, m.y.d.2)
ANDE.3 <- ANDE.fun(m.m.d.3, m.y.d.3)
ANDE.4 <- ANDE.fun(m.m.d.4, m.y.d.4)
ANDE.5 <- ANDE.fun(m.m.d.5, m.y.d.5)

c.1 <- apply(ANDE.1, 2, mean)
c.2 <- apply(ANDE.2, 2, mean)
c.3 <- apply(ANDE.3, 2, mean)
c.4 <- apply(ANDE.4, 2, mean)
c.5 <- apply(ANDE.5, 2, mean)

s.1 <- apply(ANDE.1, 2, sd)
s.2 <- apply(ANDE.2, 2, sd)
s.3 <- apply(ANDE.3, 2, sd)
s.4 <- apply(ANDE.4, 2, sd)
s.5 <- apply(ANDE.5, 2, sd)

coefs <- cbind(c.1, c.2, c.3, c.4, c.5)
c.combined <- apply(coefs, 1, mean)

var.within <- cbind(s.1^2, s.2^2, s.3^2, s.4^2, s.5^2)
var.within <- apply(var.within, 1, mean)

m <- 5
var.between <- (coefs - c.combined)^2
var.between <-  apply(var.between, 1, sum)
var.between <-  (m-1)^-1 * var.between

s.combined <- sqrt(var.within + (1 + m^-1) * var.between )

combined.res <- cbind(c.combined, s.combined)
combined.res <- round(combined.res, 2)

av.ANDE.off <- combined.res[3,1]
av.ANDE.se.off <- combined.res[3,2]
